Tarea 4: Bayesian Inference

CC6104: Statistical Thinking

Integrantes :

Tamara Marciel
Vicente Pinochet

Cuerpo Docente:

Profesor: Felipe Bravo M.
Auxiliar: Martín Paredes, Benjamín Farías
Ayudantes: Scarleth Betancurt, Emilio Díaz, Kevin Iturra, Renzo Zanca

Índice:

Objetivo
Instrucciones
Referencias
Elaboración de Código

Objetivo

Bienvenides a la cuarta tarea del curso Statistical Thinking. Esta tarea tiene como objetivo evaluar los contenidos teóricos de la primera parte del curso, los cuales se enfocan principalmente en introducirlos en la estadística bayesiana. Si aún no han visto las clases, se recomienda visitar los enlaces de las referencias.

La tarea consta de una parte práctica con el fin de introducirlos a la programación en R enfocada en el análisis estadístico de datos.

Instrucciones:

La tarea se realiza en grupos de máximo 2 personas. Pero no existe problema si usted desea hacerla de forma individual.
La entrega es a través de u-cursos a más tardar el día estipulado en la misma plataforma. A las tareas atrasadas se les descontará un punto por día.
El formato de entrega es este mismo Rmarkdown y un html con la tarea desarrollada. Por favor compruebe que todas las celdas han sido ejecutadas en el archivo html.
Al momento de la revisión tu código será ejecutado. Por favor verifica que tu entrega no tenga errores de compilación.
No serán revisadas tareas desarrolladas en Python.
Está PROHIBIDO la copia o compartir las respuestas entre integrantes de diferentes grupos.
Pueden realizar consultas de la tarea a través de U-cursos y/o del canal de Discord del curso.

Referencias:

Slides de las clases:

Videos de las clases:

Introduction to Bayesian Inference: video1 video2 video3 video4 video5 video6 video7
Summarizing the Posterior: video1 video2 video3 video4

Documentación:

#Elaboración de Código

En la siguiente sección deberá resolver cada uno de los experimentos computacionales a través de la programación en R. Para esto se le aconseja que cree funciones en R, ya que le facilitará la ejecución de gran parte de lo solicitado.

Para el desarrollo preste mucha atención en los enunciados, ya que se le solicitará la implementación de métodos sin uso de funciones predefinidas. Por otro lado, Las librerías permitidas para desarrollar de la tarea 4 son las siguientes:

# Manipulación de estructuras
library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(dplyr)
library(tidyr)
library(purrr)

# Para realizar plots
library(scatterplot3d)
library(ggplot2)
library(plotly)

## 
## Adjuntando el paquete: 'plotly'
## 
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## 
## The following object is masked from 'package:stats':
## 
##     filter
## 
## The following object is masked from 'package:graphics':
## 
##     layout

# Manipulación de varios plots en una imagen.
library(gridExtra)

## 
## Adjuntando el paquete: 'gridExtra'
## 
## The following object is masked from 'package:dplyr':
## 
##     combine

# Análisis bayesiano
library("rethinking")

## Cargando paquete requerido: rstan

## Warning: package 'rstan' was built under R version 4.4.2

## Cargando paquete requerido: StanHeaders

## Warning: package 'StanHeaders' was built under R version 4.4.2

## 
## rstan version 2.32.6 (Stan version 2.32.2)
## 
## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)
## For within-chain threading using `reduce_sum()` or `map_rect()` Stan functions,
## change `threads_per_chain` option:
## rstan_options(threads_per_chain = 1)
## 
## Do not specify '-march=native' in 'LOCAL_CPPFLAGS' or a Makevars file
## 
## Adjuntando el paquete: 'rstan'
## 
## The following object is masked from 'package:tidyr':
## 
##     extract
## 
## Cargando paquete requerido: parallel
## Cargando paquete requerido: dagitty
## rethinking (Version 2.01)
## 
## Adjuntando el paquete: 'rethinking'
## 
## The following object is masked from 'package:purrr':
## 
##     map
## 
## The following object is masked from 'package:stats':
## 
##     rstudent

Si no tiene instalada la librería “rethinking”, ejecute las siguientes líneas de código para instalar la librería:

install.packages("rethinking")

## Warning: package 'rethinking' is in use and will not be installed

En caso de tener problemas al momento de instalar la librería con el código anterior, utilice las siguiente chunk:

install.packages("shape")
install.packages(c("mvtnorm","loo","coda"), repos="https://cloud.r-project.org/",dependencies=TRUE)
options(repos=c(getOption('repos'), rethinking='http://xcelab.net/R'))
install.packages('rethinking',type='source')

Pregunta 1:

Un conjunto de carteros aburridos de las mordidas de perros ha decidido realizar un catastro de mordidas recibidas por los empleados de su empresa en un periodo de 8 meses, planeando en base a estos datos realizar inferencia bayesiana. Los datos de las mordidas estas datos por el dataset no+mordidas.csv, en donde cada fila representa las mordidas recibidas por diferentes carteros y las columnas señalan si fue mordido o no el cartero en los meses de estudio (si fue mordido es señalado con un 1, de lo contrario es señalado con un 0). Cabe señalar que un cartero no puede ser mordido mas de una vez al mes, ya que el damnificado recibe licencia de todo un mes tras la mordida, reincorporándose el siguiente mes al trabajo.

dataMordidas <- read.csv("no_mordidas.csv", header = TRUE)

[] Uno de los carteros le comenta que, según él, lo han mordido 6 de cada 10 veces, y le recomienda usar esta información como antecedente para su inferencia. Para ello considere como prior la distribución gausseana de media 0.6 y desviación estándar 0.05. Programe el método grid approximation para 1, 2, 4 y 8 meses. ¿Como van cambiando las curvas posterior?

cols <- c("bites_month_1", "bites_month_2", "bites_month_3",
          "bites_month_4", "bites_month_5", "bites_month_6",
          "bites_month_7", "bites_month_8")

grid_len <- 100
p_grid <- seq( from=0, to=1, length.out= grid_len)
prior <- dnorm(p_grid, 0.6, 0.05)

posteriors <- c()
counted.months <- c()
p_grids <- rep(p_grid, 4)

for (n in c(1, 2, 4, 8)) {
  x <- sum(dataMordidas[, cols[1:n]])# total de mordidas
  
  likelihood <- dbinom(x, n * nrow(dataMordidas), p_grid)
  unstd.posterior <- likelihood *prior
  posterior <- unstd.posterior / sum(unstd.posterior)
  
  posteriors <- c(posteriors, posterior)
  counted.months <- c(counted.months, rep(n, grid_len))
}

post.df <- data.frame("grid"=p_grids, "posterior"=posteriors, "months"=counted.months)

ggplot(post.df) +
  geom_line(aes(x=grid, y=posterior)) +
  facet_grid(months ~ .)

Respuesta:

Las curvas van adquiriendo una mayor altura y ocurriendo en un rango menor de valores para el eje x, esto supone que los datos que son usados para el modelo se van concentrando cada vez más en un rango más acotado de valores.

[] Vuelva a aplicar grid approximation como el parte anterior, pero esta vez con prior uniforme y compare los resultados. ¿Qué puede decir sobre la elección de los priors?

prior <- rep(1, grid_len)
posteriors <- c()

for (n in c(1, 2, 4, 8)) {
  x <- sum(dataMordidas[, cols[1:n]])
  
  likelihood <- dbinom(x, n * nrow(dataMordidas), p_grid)
  unstd.posterior <- likelihood *prior
  posterior <- unstd.posterior / sum(unstd.posterior)
  
  posteriors <- c(posteriors, posterior)
}

post.df$posterior2 <- posteriors

ggplot(post.df) +
  geom_line(aes(x=grid, y=posterior)) +
  geom_line(aes(x=grid, y=posterior2), color="red") +
  facet_grid(months ~ .)

Respuesta:

Las distribuciones para el prior uniforme están más a la izquierda que las de la prior gaussiana, esto ocurre ya que el prior uniforme asigna la misma probabilidad a todos lo datos, en cambio con el prior gaussiano se hace más probable el hecho de ser mordido por un perro y es por esto que se encuentra más a la derecha.

[] Repita el mismo analisis anterior pero utilizando el metodo de Laplace (no necesita programar el metodo, puede utilizar la libreria “rethinking”). ¿Cómo se comparan con los resultados anteriores?

 for (n in c(1, 2, 4, 8)) {
  res <- quap(
    alist(
      W ~  dbinom (W+L, p) ,  # binomial likelihood
      p ~  dunif(0,1)# uniform pior
    ),
    data = list(W=sum(dataMordidas[, cols[1:n]]),L=n*250-sum(dataMordidas[, cols[1:n]])) )
  
  df <- precis(res)
  curve(dnorm(x, df$mean, df$sd), lty=2, add=FALSE)
}

Respuesta:

En este caso para todos los rangos de meses las distribuciones se ven centradas en un valor de 0.4 y al utilizar más datos para el modelo se van viendo menos casos lejos del centro, cada vez con “colas” más pequeñas, esto sucede ya que el metodo de laplace solo ve la punta que distribuye como normal y por tanto se tendrán valores más acotados, es decir la curva ocurrirá en rangos más “delgados”.

[] Grafique la densidad de la posterior y encuentre la proporción de los siguientes defined boundaries:
\((0, 0.35)\)
\((0.35, 0.45)\)
\((0.45, 1)\)

¿Cómo puede interpretar los resultados?

grid_len <- 1000
p_grid <- seq( from=0, to=1, length.out= grid_len)

prior <- rep(1, grid_len)

likelihood <- dbinom(6, 10 , p_grid)

unstd.posterior <- likelihood *prior

posterior <- unstd.posterior/sum(unstd.posterior)

samples <- sample(p_grid, prob=posterior, size=1e4, replace=TRUE)

dens(samples)

Intervalos:

# fracción de valores <0.35
sum(posterior [p_grid < 0.35])

## [1] 0.05002864

# fracción de valores entre 0.35 y 0.45
sum(posterior [0.35 < p_grid & p_grid < 0.45])

## [1] 0.1236835

# fracción de valores >0.45
sum(posterior [ 0.45 < p_grid])

## [1] 0.8262878

Respuesta:

Cerca de un 5% del área de probabilidad del posterior está bajo un valor de 0.35.
Cerca de un 12,4% del área de probabilidad del posterior está entre un valor de 0.35 y un valor de 0.45.
Cerca de un 82,6% del área de probabilidad del posterior está sobre un valor de 0.45.

Lo cual sugiere que la concentración de los resultados está ligeramente inclinada hacia la derecha, tal y como se ve en el gráfico de la distribución de la posterior.

-[] Calcule un intervalo de credibilidad al \(50%\), \(75%\) y \(95%\). ¿Como se puede interpretar los resultados? ¿Cual podría ser un problema al usar intervalos de credibilidad?

Intervalo \(50\%\):

quantile(samples, c(0.25, 0.75))

##       25%       75% 
## 0.4884885 0.6826827

PI(samples, 0.5)

##       25%       75% 
## 0.4884885 0.6826827

Intervalo \(75\%\)

quantile(samples, c(0.125, 0.875))

##     12.5%     87.5% 
## 0.4184184 0.7427427

PI(samples, 0.75)

##       12%       88% 
## 0.4184184 0.7427427

Intervalo \(95\%\)

quantile(samples, c(0.025, 0.975))

##      2.5%     97.5% 
## 0.3103103 0.8278529

PI(samples, 0.95)

##        3%       98% 
## 0.3103103 0.8278529

Respuesta:

Hay un 50% de probabilidad de obtener un valor entre 0.49 y 0.68.
Hay un 75% de probabilidad de obtener un valor entre 0.42 y 0.74.
Hay un 95% de probabilidad de obtener un valor entre 0.3 y 0.83.

Es un rango de valores que creemos que el parámetro va a tomar con una cierta probabilidad.

Un problema puede surgir si las distribuciones obtenidas del posterior se comportan de una manera muy asimetrica, tal que las colas arrojen resultados engañosos y no congruentes a los intervalos requeridos.

Genere los intervalos HDPI para \(50%\), \(75%\) y \(95%\), compárelos con los intervalos de credibilidad de la parte anterior. ¿En que se diferencian los intervalos de credibilidad con los HDPI?

Intervalo \(50\%\):

HPDI(samples, 0.5)

##      |0.5      0.5| 
## 0.4874875 0.6816817

Intervalo \(75\%\)

HPDI(samples, 0.75)

##     |0.75     0.75| 
## 0.4284284 0.7507508

Intervalo \(95\%\)

HPDI(samples, 0.95)

##     |0.95     0.95| 
## 0.3133133 0.8298298

Respuesta:

Dio como resultado valores bastante similares a los resultados de la parte anterior, esto era esperable ya que la distribución es bastante simétrica y por tanto las “colas” generadas en ambos serían bastante similares ya que la mayor concentración está en el centro.

los intervalos de HPDI buscan los intervalos más delgados y no necesariamento los que contengan la misma proporción de valores en las “colas”, además son más costosos.

Utilizando la posterior obtenida en la parte anterior, simule 10.000 réplicas de registros de mordidas. Compare la distribución del ratio de carteros mordidos predichos a partir de las réplicas con el ratio real de los datos. ¿El modelo se ajusta bien a los datos? Es decir, ¿la distribución de las predicciones incluye la observación real como resultado central y probable?

# Número total de carteros y datos observados
num_carteros <- nrow(dataMordidas)
num_meses <- 8
total_mordidas <- sum(dataMordidas[, cols])  # total de mordidas

# Ratio real de mordidas
ratio.real <- total_mordidas 
samples.sim <-  rbinom(1e4, num_carteros*num_meses, prob= samples)

ggplot() +
  geom_density(aes(samples.sim)) +
  geom_vline(aes(xintercept=ratio.real), color="red")

Respuesta:

Este modelo no se ajusta bien a los datos, ya que la distribución de las predicciones no incluye la observación real como resultado central y probable, esto sucede ya que se tomó como hecho más probable el ser mordido, siendo la realidad que la frecuencia de carteros mordidos es mucho menor a lo asumido.

Escoja cualquier par de meses concecutivos y obtenga la posterior de que una persona que fue mordida en el primer mes, sea mordida nuevamente en el segundo mes. Para esto, se recomienda comenzar buscando los carteros que fueron mordidos el primer mes y en base a estos generar una búsqueda indexada para obtener el número solicitado. Hecho esto simule ese número carteros mordidos 10.000 veces. De los resultados obtenidos, compare el recuento de carteros mordidos con el recuento real. ¿Cómo se ve el modelo desde este punto de vista?

# Se escogen el mes 5 y 6
primer_mes <- dataMordidas$bites_month_5
segundo_mes <- dataMordidas$bites_month_6

# Carteros que fueron mordidos en el primer mes
mordidos_primer_mes <- which(primer_mes == 1)
cantidad_mordidos_primer_mes <- length(mordidos_primer_mes)

# Carteros mordidos en el segundo mes dado que fueron mordidos en el primero
mordidos_condicionados <- sum(segundo_mes[mordidos_primer_mes])

# Calculando la posterior con prior uniforme
grid_len <- 1000
p_grid <- seq( from=0, to=1, length.out= grid_len)

prior <- rep(1, grid_len)

likelihood <- dbinom(mordidos_condicionados, cantidad_mordidos_primer_mes , p_grid)

unstd.posterior <- likelihood *prior

posterior <- unstd.posterior/sum(unstd.posterior)

samples <- sample(p_grid, prob=posterior, size=1e4, replace=TRUE)

# Simulación de 10,000 réplicas

simulaciones <- rbinom(1e4, size = cantidad_mordidos_primer_mes, prob = samples)

ggplot() +
  geom_density(aes(simulaciones)) +
  geom_vline(aes(xintercept=mordidos_condicionados), color="red")

Respuesta:

Se puede ver una distribución bastante simétrica centrada en el valor del recuento real, por lo que se puede decir que son simulaciones bastante apegadas a la realidad de los datos del modelo.

Pregunta 2:

En esta pregunta se trabajara con el dataset “notas.csv” el cual contiene las notas históricas de un curso desconocido. Suponga que los datos vienen de una distribución \(\mathcal{N}(\mu,\sigma^2)\), el objetivo de la pregunta es estudiar el comportamiento de los datos y los posibles valores de \(\mu,\sigma\) mediante técnicas de inferencia Bayesiana.

Usted sabe un dato extra sobre la información, \(\sigma \in [0.5,1.5]\) y, además, tiene una fuerte creencia a que es mas probable encontrar la desviación estándar real entre \([0.5,1]\) que en \((1,1.5]\).

Modifique el siguiente codigo que permite realizar una grid approximation para \(2\) parámetros. Proponga priors para \(\mu\) y \(\sigma\), justifique su elección.

# Leer información
data_notas <- read.csv("notas.csv")
head(data_notas)

# Función para crear likelihood dado mu y sigma
grid_function <- function(mu, sigma) {
  # Log-likelihood basado en distribución normal
  sum(dnorm(data_notas$Notas, mean = mu, sd = sigma, log = TRUE))
}

# Valores de la grilla
grid_mu <- seq(mean(data_notas$Notas, na.rm = TRUE) - 3, 
               mean(data_notas$Notas, na.rm = TRUE) + 3, length.out = 100)
grid_sigma <- seq(0.5, 1.5, length.out = 100)

# Crear grilla 2D
data_grid <- expand_grid(grid_mu, grid_sigma)

# Calcular likelihood para cada par (mu, sigma)
data_grid <- data_grid %>%
  mutate(likelihood = map2_dbl(grid_mu, grid_sigma, grid_function))

# Priors
prior_mu <- rep(1 / length(grid_mu), length(grid_mu)) # Uniforme para mu
prior_sigma <- ifelse(grid_sigma <= 1, 0.75, 0.25)    # Ponderación más alta para [0.5, 1]
prior_sigma <- prior_sigma / sum(prior_sigma)         # Normalizar prior_sigma

# Crear columnas separadas para priors en data_grid
data_grid <- data_grid %>%
  mutate(prior_mu = prior_mu[match(grid_mu, unique(grid_mu))],
         prior_sigma = prior_sigma[match(grid_sigma, unique(grid_sigma))],
         prior = prior_mu * prior_sigma)

# Calcular posterior no estandarizado
data_grid <- data_grid %>%
  mutate(unstd_posterior = likelihood + log(prior), # Log-transformación para estabilidad numérica
         posterior = exp(unstd_posterior - max(unstd_posterior))) %>% # Ajuste para evitar underflow
  mutate(posterior = posterior / sum(posterior)) # Normalización del posterior

# Mostrar las primeras filas del dataframe resultante
head(data_grid)

Respuesta:

Prior para \(\mu\) (media):
Dado que no tenemos información previa sobre la ubicación de la media, se asume una distribución uniforme para \(\mu\). Esto refleja una incertidumbre inicial equitativa sobre los posibles valores de la media.
Prior para \(\sigma\) (desviación estándar):
Sabemos que \(\sigma \in [0.5, 1.5]\), con una mayor probabilidad en \([0.5, 1]\) que en \((1, 1.5]\).
Para reflejar esta creencia, asignamos un prior con un peso de 0.75 para valores en \([0.5, 1]\) y 0.25 para valores en \((1, 1.5]\). Este prior sesgado refleja nuestra expectativa sobre \(\sigma\).

Tras haber ejecutado el código de la parte anterior ejecute el siguiente, ¿Que puede decir de los valores de la distribución? Se recomienda hacer modificaciones en el código para realizar un mejor análisis y estudio.

# Punto de referencia
# Se recomienda cambiar estos valores por unos adecuados que le permitan estudiar
# Los valores de la distribución de mejor manera
valor_x <- 1
valor_y <- 1

# Grafico

punto_comparacion <- tibble(x = valor_x, y = valor_y)

plt <- data_grid %>%
  ggplot(aes(x = grid_mu, y = grid_sigma)) +
  geom_raster(aes(fill = posterior),
    interpolate = T
  )+ 
  geom_point(x = valor_x, y = valor_y, size = 1.3,color="white")+
  geom_label(
    data = punto_comparacion, aes(x, y),
    label = "Punto Comparación",
    fill = "green",
    color = "black",
    nudge_y = 0, # Este parametro desplaza la caja por el eje y
    nudge_x = 1 # Este parametro desplaza la caja por el eje x
  )+
  scale_fill_viridis_c() +
  labs(
    title = "Posterior para Mean y Standard Deviation",
    x = expression(mu ["Mean"]),
    y = expression(sigma ["Standar Deviation"])
  ) +
  theme(panel.grid = element_blank())

plt

Respuesta:

Observaciones sobre la distribución posterior

Zona de máxima densidad: La mayor probabilidad para los parámetros \(\mu\) (media) y \(\sigma\) (desviación estándar) se concentra en un área definida. Esto indica que la estimación conjunta de los parámetros es consistente con valores específicos, sugiriendo un comportamiento estable en la media y la desviación estándar.
Forma de la distribución: La alta concentración de probabilidad en una región angosta sugiere baja incertidumbre sobre la estimación de \(\mu\) dentro del rango observado. La distribución está sesgada hacia un intervalo reducido para \(\sigma\), con la mayor densidad posterior centrada en valores cercanos a 0.75-1, lo que coincide con la creencia previa de que \(\sigma\) se encuentra en el rango \([0.5, 1]\).

A continuación se presenta un código que permite realizara la distribución dada por sampling from grid approximation ¿Para que sirve este proceso? ¿Que puede deducir del gráfico?

# Codificamos los datos
x <- 1:length(data_grid$posterior)

# Sampleamos los indices
posterior_samples_aux <- sample(x,size = 1e4, replace = T, prob = data_grid$posterior)

# Obtenemos los verdaderos valores de la sampling distribution
posterior_samples <- data_grid[posterior_samples_aux,]

# Obtenemos solos los valores relevantes para la densidad
df <- data.frame(posterior_samples$grid_mu,posterior_samples$grid_sigma)

# Realizamos las densidades
dens(df)

Respuesta:

El proceso de sampling from grid approximation sirve para estudiar cómo se distribuyen los valores de los parámetros \(\mu\) y \(\sigma\) en la posterior. En el gráfico, se pueden observar las regiones de mayor densidad de puntos, lo cual sugiere que los valores de \(\mu\) y \(\sigma\) más consistentes con los datos se concentran en esas áreas. El gráfico muestra una mayor densidad en la parte baja del rango de \(\sigma\) (entre 0.55 y 0.7), esto confirmaría la prior que propusimos, que da mayor probabilidad a esa región. Además, la forma y dispersión de los puntos indican la certeza o incertidumbre en las estimaciones de los parámetros.

En cuanto a la distribución de la media (\(\mu\)), podemos observar que los valores más probables se encuentran en la región central del gráfico, donde se concentra la mayor densidad de puntos. Esto sugiere que la media estimada se encuentra dentro del rango 5.9 y 6.2, el cual indica la consistencia de los datos con los valores propuestos en la prior. Los puntos cercanos al centro de la grilla reflejan una mayor certeza sobre el valor de \(\mu\), mientras que aquellos más alejados de este centro reflejan mayor incertidumbre.

Pregunta 3: Regresión Lineal Bayesiana

El objetivo de esta pregunta es introducirlo en la aplicación de una regresión bayesiana. Con esto, buscaremos entender como calcular una regresión bayesiana en base al “motor” aproximación de Laplace, revisando como se comporta la credibilidad de sus predicciones y como la regresión lineal puede llegar a mutar a aplicar una transformación en el vector \(x\). Para responder esta pregunta centre su desarrollo solo en las clases y las funciones que nos brinda la librería rethinking.

Para esta pregunta usaremos el dataset Pokemon.csv, que contiene las características de más de 700 personajes.

poke.df <- read.csv("Pokemon.csv")
head(poke.df)

Defina un modelo de regresión bayesiana para modelar Sp..Atk en función de Attack, definiendo sus propios priors. Comente cada una de sus asunciones y señale a través de ecuaciones el modelo.
Ajuste el modelo lineal utilizando el método de Laplace approximation. Estudie a través del comando precis los resultados obtenidos y coméntelos.
Gráfique el MAP de regresión lineal obtenida, añadiendo al gráfico el intervalo del \(95\%\) que se tiene sobre la media y los valores predecidos de la altura. Comente los resultados obtenidos y señale si su modelo logra ser un buen predictor de los valores estudiados.

Respuesta:

modelo <- alist(
  Sp..Atk ~ dnorm(mu, sigma),
  mu <- alpha + beta * Attack,
  alpha ~ dnorm(100, 20),   # Prior para alpha
  beta ~ dnorm(0, 10),      # Prior para beta
  sigma ~ dexp(1)           # Prior para sigma
)

El modelo propuesto busca modelar la relación entre Sp..Atk (ataque especial) y Attack (ataque físico). Las ecuaciones y los priors asumidos son los siguientes:

Modelo matemático:

\[ Sp..Atk_i \sim \mathcal{N}(\mu_i, \sigma) \quad \text{donde} \quad \mu_i = \alpha + \beta \cdot \text{Attack}_i \]

Priors definidos: - \(\alpha \sim \mathcal{N}(100, 20)\): Refleja que esperamos un valor promedio de Sp..Atk alrededor de 100 con una incertidumbre razonable. - \(\beta \sim \mathcal{N}(0, 10)\): Asume inicialmente que no hay una relación fuerte entre Attack y Sp..Atk. - \(\sigma \sim \text{Exponencial}(1)\): Prior para la desviación estándar de los residuos, asegurando positividad y favoreciendo valores pequeños.

fit <- map(
  alist(
    Sp..Atk ~ dnorm(mu, sigma),
    mu <- alpha + beta * Attack,
    alpha ~ dnorm(100, 20),
    beta ~ dnorm(0, 10),
    sigma ~ dexp(1)
  ),
  data = poke.df
)
 
precis(fit)

El intercepto (\(\alpha\)) indica que, para un valor de Attack cercano a 0, el Sp..Atk predicho es 42.41. La pendiente (\(\beta\)) es 0.39, lo que implica que por cada aumento de una unidad en Attack, el Sp..Atk promedio aumenta en 0.39 unidades. El alto valor de \(\sigma\) (29.50) refleja una gran variabilidad no explicada por el modelo.

plot_data <- data.frame(Attack = seq(min(poke.df$Attack), max(poke.df$Attack), length.out = 100))
plot_data$mu <- coef(fit)["alpha"] + coef(fit)["beta"] * plot_data$Attack
ggplot() +
  geom_point(aes(x = poke.df$Attack, y = poke.df$Sp..Atk), color = "red") +
  geom_line(aes(x = plot_data$Attack, y = plot_data$mu), color = "orange") +
  geom_ribbon(aes(x = plot_data$Attack, ymin = plot_data$mu - 1.96 * coef(fit)["sigma"], 
                  ymax = plot_data$mu + 1.96 * coef(fit)["sigma"]), alpha = 0.3, fill = "orange")

Comentarios sobre el gráfico

Relación positiva moderada:
- La línea naranja representa el modelo de regresión ajustado. Se observa una tendencia positiva entre Attack y Sp..Atk, lo que confirma que al aumentar el valor de Attack, también aumenta el valor promedio de Sp..Atk.
Alta dispersión en los datos:
- Los puntos rojos (datos reales) están altamente dispersos alrededor de la línea de regresión. Esto indica que el modelo no captura toda la variabilidad en los datos. La dispersión sugiere que existen otros factores (además de Attack) que influyen significativamente en Sp..Atk.
Intervalos de credibilidad amplios:
- Las áreas sombreadas en naranja representan los intervalos de credibilidad al 95%. La amplitud considerable de estos intervalos refuerza la idea de alta incertidumbre en las predicciones del modelo.
Valores atípicos:
- Hay puntos alejados significativamente de la línea ajustada, lo que podría corresponder a Pokémon con características muy específicas (por ejemplo, Pokémon legendarios o de tipos especializados).
1. Limitaciones del modelo lineal:
  - Aunque el modelo lineal captura la tendencia general, parece insuficiente para describir completamente la relación entre Attack y Sp..Atk. Un modelo más complejo podría mejorar la capacidad predictiva.
  El modelo muestra que existe una relación positiva, pero no captura toda la variabilidad ni explica los datos de manera óptima debido a la alta dispersión y otros factores omitidos.

A work by CC6104